####################################
# Covariates RDD
####################################

rm(list=ls())

library(Hmisc)
library(ggplot2)
library(stargazer)
library(foreign)
library(rdrobust)
library(rdd)
library(readstata13)
library(haven)

################
# Prepare data 
################

# read data
d = read_dta("~/Dropbox/Anti-american attitudes/01_data/clean_data/trump_election_data_clean_final.dta")   
names(d)

# bandwidths
h_female = rdbwselect(d$female,d$score)$bws[1]
h_female

h_age = rdbwselect(d$age,d$score)$bws[1]
h_age

h_education = rdbwselect(d$education,d$score)$bws[1]
h_education

h_cluster = rdbwselect(d$cluster2,d$score)$bws[1]
h_cluster

h_left = rdbwselect(d$left,d$score)$bws[1]
h_left

######################
# Female
######################

kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$female ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$female ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$female ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$female))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_female))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)
results$pv[results$color=="red"]
mean(results$pv)

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(legend.position="none")

######################
# Age
######################

kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$age ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$age ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$age ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$age))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_age))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)
results$pv[results$color=="red"]
mean(results$pv)

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_y_continuous(breaks=c(-4,-3.5,-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4)) + coord_cartesian(ylim=c(-4,4)) + geom_hline(yintercept=0) + xlab("Bandwidth") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(legend.position="none")

######################
# Education
######################

kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$education ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$education ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$education ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$education))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_education))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)
results
results$pv[results$color=="red"]
mean(results$pv)

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_y_continuous(breaks=c(-0.5,0,0.5)) + coord_cartesian(ylim=c(-0.5,0.5)) + geom_hline(yintercept=0) + xlab("Bandwidth") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(legend.position="none")

######################
# Cluster
######################

kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$cluster2 ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$cluster2 ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$cluster2 ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$cluster2))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_cluster))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)
results$pv[results$color=="red"]
mean(results$pv)

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_y_continuous(breaks=c(-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10)) + coord_cartesian(ylim=c(-10,10)) + geom_hline(yintercept=0) + xlab("Bandwidth") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(legend.position="none")

######################
# Left
######################

kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$left ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$left ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$left ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$left))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_left))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)
results$pv[results$color=="red"]
mean(results$pv)

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(legend.position="none")

